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Exact Wavelets on the Ball 
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Abstract — We develop an exact wavelet transform on the three- 
dimensional ball (i.e. on the solid sphere), which we name 
the flaglet transform. For this purpose we first construct an 
exact transform on the radial half-line using damped Laguerre 
polynomials and develop a corresponding quadrature rule. Com- 
bined with the spherical harmonic transform, this approach 
leads to a sampling theorem on the ball and a novel three- 
dimensional decomposition which we call the Fourier-Laguerre 
transform. We relate this new transform to the well-known 
Fourier-Bessel decomposition and show that band-limitedness in 
the Fourier-Laguerre basis is a sufficient condition to compute 
the Fourier-Bessel decomposition exactly. We then construct the 
flaglet transform on the ball through a harmonic tiling, which is 
exact thanks to the exactness of the Fourier-Laguerre transform 
(from which the name flaglets is coined). The corresponding 
wavelet kernels are well localised in real and Fourier-Laguerre 
spaces and their angular aperture is invariant under radial 
translation. We introduce a multiresolution algorithm to perform 
the flaglet transform rapidly, while capturing afl information at 
each wavelet scale in the minimal number of samples on the 
baU. Our implementation of these new tools achieves floating- 
point precision and is made publicly available. We perform 
numerical experiments demonstrating the speed and accuracy 
of these libraries and illustrate their capabilities on a simple 
denoising example. 

Index Terms — Harmonic analysis, wavelets, ball. 



L Introduction 

A common problem in data analysis is the extraction 
of non-trivial patterns and structures of interest from 
signals. This problem can be addressed by projecting the data 
onto an appropriate basis. Whereas Fourier analysis focuses 
on oscillatory features, wavelets extract the contributions of 
scale-dependent features in both real and frequency space 
simultaneously. Initially defined in Euclidean space, wavelets 
have been extended to various manifolds and are now widely 
used in numerous disciplines. In particular, spherical wavelets 
P|-|TQ| have been extremely successful at analysing data on 
the sphere and have now become a standard tool in geophysics 
{e.g. |[TT|-|[T6|) and astrophysics {e.g. |T7|-|'32|). Naturally, 
data may also be defined on the three-dimensional ball when 
radial information (such as depth, redshift or distance, for 
example) is associated with each spherical map. 

First approaches to perform wavelet-type transforms on the 
ball were developed by p3| , p4| in the continuous setting 
only, which thus cannot be used for exact reconstruction 
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in practice. The spherical Haar transform p5| , p6 | was 
extended to the ball by [37J to support exact analysis and 
synthesis. However, this framework is very restrictive and 
may not necessarily lead to a stable continuous basis |38|- 
[40 1 . The first wavelet transform on the ball to tackle both 
the continuous and discrete settings was developed in the 
influential work of Lanusse et al. |41 1. This wavelet transform 
is based on an isotropic undecimated wavelet construction, 
built on the Fourier-Bessel transform. Since these wavelets are 
isotropic, their angular aperture depends on the distance to the 
origin. Although the wavelet transform on the ball is exact in 
Fourier-Bessel space, wavelet coefficients must be recovered 
on the ball from their Fourier-Bessel coefficients (in order to 
extract spatially localised information). However, there exists 
no exact quadrature formula for the spherical Bessel transform 
(the radial part of the Fourier-Bessel transform) |42|, and 
thus no way to perform the Fourier-Bessel transform exactly. 
Consequently, the undecimated wavelet transform on the ball is 
not theoretically exact when wavelet coefficients are recovered 
on the ball. Nevertheless, the isotropic undecimated wavelet 
transform does achieve good numerical accuracy, which may 
be sufficient for many applications [^Wavelets on the ball have 
also been discussed in geophysics by p3| , |T4| , who espoused 
a philosophy of separability in the three Cartesian coordinates 
of a ball-to-"cubed-sphere-baH" mapping, although in |T3| 
examples are shown where wavelet transforms have been 
performed on each spherical shell only but not in the radial 
direction. In ongoing work, these same authors have extended 
their approach to the ball, where the wavelet transform in 
the radial direction is tailored to seismological applications 
by honouring certain major discontinuities in the seismic 
wavespeed profile of the Earth p31 , p6| . At present, to the 
best of our knowledge, there does not exist an exact wavelet 
transform of a band-limited signal defined on the ball. 

One reason there is no exact wavelet transform on the ball 
is due to the absence of an exact harmonic transform. We 
resolve this issue by deriving an exact spherical Laguerre 
transform on the radial half-line, leading to a new Fourier- 
Laguerre transform on the ball which is theoretically exact. 
Furthermore, this gives rise to a sampling theorem on the ball, 
where all information of a band-limited function is captured in 
a finite number of samples. With an exact harmonic transform 
on the ball in hand, we construct exact wavelets through a 
harmonic tiling, which we call flaglets (since they are built 
on the Fourier-LAGuerre transform). Each wavelet kernel is 
localised in real and Fourier-Laguerre spaces, and probes a 
characteristic angular scale which is invariant under radial 
translation. Flaglets allow one to probe three-dimensional 

^The accuracy of the Fourier-Bessel transform, and thus the isotropic 
undecimated wavelet transform on the ball, may be improved by numerical 
iteration, although this can prove problematic for certain applications. 
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spherical data in position and scale simultaneously. Moreover, 
their exactness properties guarantee that the flaglet transform 
captures and preserves all the information contained in a band- 
limited signal. 

The remainder of this article is organised as follows. In 
Section [ll| we define the spherical Laguerre transform on the 
radial half-line and the Fourier-Laguerre transform on the ball. 



In Section III we construct the exact flaglet transform on the 



ball. In Section IV we present a multiresolution algorithm 
to compute the flaglet transform and evaluate our algorithms 
numerically. A simple denoising example is presented in 
Section [V| Concluding remarks are made in Section VI 



II. Harmonic Analysis on the Ball 

The aim of this section is to construct a novel three- 
dimensional transform which is appropriate for spherical co- 
ordinates and admits an exact quadrature formula. For this 
purpose, we first set out a radial one-dimensional transform 
inspired by the Laguerre polynomials and we derive a natural 
sampling scheme and quadrature rule on the radial half- 
line. We relate this novel spherical Laguerre transform to the 
spherical Bessel transform and show that the latter can be 
evaluated exactly if the signal is band-limited in the spherical 
Laguerre basis. We combine the spherical Laguerre transform 
with the spherical harmonics to form the Fourier-Laguerre 
transform on the ball, yielding a novel sampling theorem and 
an exact harmonic transform^ 

A. The spherical Laguerre transform 

The Laguerre polynomials, solutions to the Laguerre differ- 
ential equation | [43| , | [44| , are well known for their various ap- 
plications in engineering and physics, notably in the quantum- 
mechanical treatment of the hydrogen atom ||45l, as well as 
in modern optics | [46| , | [47| . They form a natural orthogonal 
basis on the interval [0, oo) {i.e. non-negative reals R+) with 
respect to an exponential weight function. In this work, since 
we use this expansion along the radial half-line, we define the 
spherical Laguerre basis function Kp{r) with r G as 



Kp{r) 



r/27 



(p + 2)! 



(1) 



(2) 

where Lp ^ is the p-th generalised Laguerre polynomial of 
order two, defined as 

and r G is a scale factor that adds a scaling flexibility 
and shall be defined at the end of this section. The basis 

^A harmonic transform is typically associated with basis functions which 
are eigenvalues of the Laplacian operator (e.g. the Fourier transform). In 
this paper our basis functions on the radial half-line (and thus on the ball) 
are not solutions of the Laplacian, hence harmonic analysis on the ball 
is interpreted in a broader sense. Nonetheless, these basis functions form 
orthonormal transforms and define valid dual spaces. We define band-limited 
signals to have bounded support in the transform space of these orthogonal 
basis functions. 



functions Kp are orthonormal on R+ with respect to a radial 
inner product: 



{Kp\K,) = [ drr^Kpir)K;ir) = 5p,. 



(3) 



Note that the complex conjugate * is facultative since we 
use real basis functions. Any square-integrable real signal 
/ G L^(R+) may be expanded as 



f{r) = ^f,K,{r), 

p=0 



(4) 



for natural p G N, where fp is the projection of / onto the 
p-th basis function: 

U = {f\Kp) = f drr'f{r)K;{r). (5) 

The decomposition follows by the orthogonality and com- 
pleteness of the spherical Laguerre basis functions: orthonor- 
mality is given by Eqn. ([3]), while the completeness relation 
is obtained by applying the Gram-Schmidt orthogonalisation 
process to the basis functions and exploiting the completeness 
of polynomials on L^(R+, r^e~^dr). 

When it comes to calculating the transform, one must 
evaluate the integral of Eqn. ^ numerically. We consider 
functions / band-limited at P in the spherical Laguerre basis, 
such that /p = 0, > P. It is straightforward to show that 
if / is band-limited, then both e^^'^^Kp{r) and e^/^'^/(r) are 
polynomials of maximum degree P— 1. In this case, Eqn. ([5]) is 
the integral of a polynomial of order 2P— 2 on with weight 
function r^e~^. Thus, applying Gaussian quadrature {e.g. | [48| , 
| [49| ) with P sampling nodes is sufficient to evaluate this 
integral exactly. The resulting quadrature formula is known 
as the Gauss-Laguerre quadrature and is commonly used to 
evaluate numerical integrals on R+. Hence, Eqn. ^ reduces 
to a weighted sum: 



p-i 



fp=^ Wif{ri)K^{ri), 



(6) 



i=0 



where G R^ is the i-th root of the P-th generalised 
Laguerre polynomial of order two, and 



Wi 



(P + l)[4') (r,)]2 



(7) 



is the corresponding weight. Any P-band-limited function 
/ can be decomposed and reconstructed exactly using the 
spherical Laguerre transform. All information content of the 
function is captured in P samples located in the interval 
[0,rp_i] where rp_i is the largest root of the sampling. 
Since rp_i increases with P, one may wish to rescale the 
sampling so that the spherical Laguerre transform contains 
samples in any interval of interest [0, i?], with R G R+, while 
the underlying continuous function is nevertheless defined on 
R+. The scale factor r is then chosen such that r = R/rp-i. 
Figure [T] shows the resulting spherical Laguerre sampling 
constructed on r G [0, 1] {i.e. rescaled with r) for increasing 
band-limit P. Figure |2] shows the first six basis functions 
constructed on r G [0, 1] and the sampling nodes used to 
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Fig. L Spherical Laguerre sampling scheme on r G [0, 1] for increasing 
band-limit P. If a function / is P-band-limited then / and the basis functions 
need only be evaluated on P points for the spherical Laguerre transform to 
be exact. For a particular P, the associated sampling is denser near the origin 
since the quadrature is constructed on R+ with measure e~^dr. 
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(a) Basis functions Kp(r) 




(b) Zoom on the oscillatory features of Kp(r) 



obtain an exact transform|^ Note that the spherical Laguerre 
transform is a real transform that can be extended to complex 
signals by considering the real and imaginary parts separately. 

B. Relation to the spherical Bessel transform 

The spherical Bessel transform is a fundamental radial 
transform arising from the resolution of the Laplacian operator 
in spherical coordinates. It is central to the Fourier-Bessel 
transform, commonly used in cosmology p0|-|52| to analyse 
the spectral properties of galaxy surveys in three dimensions. 
In this section we derive an analytical formula to exactly 
compute the spherical Bessel transform of a function whose 
spherical Laguerre transform is band-limited. This section is 
optional to the reader interested in wavelets only. 

The spherical Bessel transform of / G L^(R+) reads 



h{k) = {f\k) 



(8) 



for k G ^ G N, and where ji{kr) is the ^-th order spherical 
Bessel function. Note again that the complex conjugate is 
facultative since the spherical Bessel functions are real. The 
reconstruction formula is given by 



dke~f,{k)j,{kr). 



(9) 



The spherical Bessel transform is thus symmetric and the 
problem is reduced to the calculation of a similar inner product 
for the decomposition and the reconstruction. However, to 

^If one preferred to consider the measure dr rather than the spherical 
measure r^dr, then the basis functions rKp{r) shown in Figure [21(c) could 
be used in place of the spherical Laguerre basis functions defineonere. 




(c) Functions rKp(r) 

Fig. 2. First six spherical Laguerre basis functions Kp(r) constructed on 
r G [0,1] and the associated sample positions (circles). A function / with 
band-limit P = 6 can be decomposed and reconstructed exactly using these 
six basis functions only. In that case, / and the basis functions are solely 
evaluated at the sampling points. Functions rKp(r) can be viewed as basis 
functions in cartesian coordinates satisfying the usual orthogonality relation 
4+ dr(rKp(r))(rK,(r)) = 6p,. 



our knowledge, there exists no method to compute such an 
integral exactly for a useful class of functions, and finding a 
quadrature formula for the spherical Bessel functions on R+ is 
a non-trivial issue. Moreover, the use of numerical integration 
methods does not always guarantee good accuracy because of 
the oscillatory nature of the spherical Bessel functions. 

To find a tractable expression to compute Eqn. ([5]), we first 
express / by its spherical Laguerre expansion, giving 



which is a finite sum if / is band-limited in spherical Laguerre 
space. In this expression jip{k) is the projection of Kp onto 
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jeikr), i.e. 



{k) = {Kp\je)= [ drr^Kp{r)f,{kr). (11) 

Consequently, the problem of computing the spherical Bessel 
decomposition of / is recast as evaluating Eq n. (TO^ , through 
the computation of the inner product of Eqn. (111]). But unlike 
the initial problem of Eqn. ([s]), jip{k) admits an analytic 
formula. Starting from the definition of Laguerre polynomials 
in Eqn. ([2]), one can show that 



pi 



(p + 2)!f 



where the satisfy the following recurrence: 



j! \P-J 



j(j+2) ^ 

The functions are the moments of ji{kr)e~^ , i.e. 

lA(^) = \i I drr^ji{kr)e~^ . (14) 



decomposition and reconstruction operations can be performed 
with a finite summation over the harmonics. This is usually 
resolved by defining an appropriate sampling theorem on the 
sphere with nodes ujj = {Oj^(j)j), associated with a quadrature 
formula. Various sampling theorems exist in the literature 
p4|-|[56|; the main features of all sampling theorems are (i) 
the number of nodes required to capture all information in 
a band-limited signal and (ii) the complexity of the related 
algorithms to compute forward and inverse spherical harmonic 
transforms. Although this work is independent from this choice 
p (provided that it leads to an exact transform), we adopt the 

^^^2(^^)5 (12) McEwen & Wiaux (hereafter MW) sampling theorem |56| 

which is equiangular and has the lowest number of samples for 
a given band-limit L, namely (L - 1)(2L - 1) + 1 2L^. The 
corresponding algorithms to compute the spherical harmonic 
transforms scale as 0{L^) and are numerically stable to band- 
limits of at least L = 4096 |36l. Further technical details are 



(13) 



From p3| we find an analytical solution for the latter integral: 
M^(fc) = v^2^- fc^ri (15) 



where k = rk is the rescaled k scale and 2^1 is the Gaussian 
hypergeometric function. Since either (j+^+l)/2 or (j+^)/2 
is a positive integer, the latter reduces to a polynomial of P 
and it is possible to compute the quantity jip{k) exactly using 
Eqn. ( p^ to ( p3] ). Consequently, the inverse spherical Bessel 
transform fi{k) may then be calculated analytically through 
Eqn. ([To]), which is computed exactly if / is band- limited in 
the spherical Laguerre basis. 

C. The spherical harmonic transform 

Whereas the spherical Laguerre transform is specifically 
designed for analysing functions on the radial half-line, the 
spherical harmonic transform is a natural choice for the 
angular part of a consistent three-dimensional analysis. For 
a function / G L^(6'^) on the two-dimensional sphere, the 
transform reads 



(16) 



£=0 m=-l 



where = (6>, (/)) G 5^ are spherical coordinates of the unit 
sphere 6*^, with colatitude G [0, tt] and longitude (j) G [0, 27r). 
Thanks to the orthogonality and completeness of the spherical 
harmonics Yim{^), the inverse transform is given by the 
following inner product on the sphere: 



fim (/l^^ri 



52 



(17) 



with surface element do; = smOdOdcj). For a function which 
is band-limited in this basis at L, i.e. firn = 0, > L, the 



provided in Section IV-B 



D. The Fourier- Laguerre transform 

We define the Fourier-Laguerre basis functions 
on = R+ X S'^ as the product of the spherical 
Laguerre basis functions and the spherical harmonics: 
^£mp(^) = Kp{r)Y£m{'^) with the 3D spherical coordinates 
r = (r, cj) G B^. The orthogonality and completeness 
of the Fourier-Laguerre basis functions follow from the 
corresponding properties of the individual basis functions, 
where the orthogonality relation is given explicitly by the 
following inner product on B^: 



{Zirnpl^i'm'p') 



B3 



dVZ,^^Z;,^,^,(r) (18) 



-'mm' ^pp' 7 

where d^r = smOdrdOdcj) is the volume element in spher- 
ical coordinates. Any three-dimensional signal / G L'^{B^) 
can be decomposed as 



P-lL-l £ 

/(^) = f^rnpZlmpir), 
p=0 £=0 m=-£ 



(19) 



with L and P the angular and radial band-limits, respectively, 
i.e. f is such that firnp = 0, > L, Vp > P. The inverse 
relation is given by the projection of / onto the basis functions: 



{r). (20) 



The Fourier-Laguerre transform may also be relate d to th e 
Fourier-Bessel transform using the results of Section H-Bj 



In practice, calculating the Fourier-Laguerre transform re- 
quires the evaluation of the integral of Eqn. ( [20] ). For this 
purpose, combining the quadrature rules on the sphere and 

^The Fourier-Laguerre and the Fourier-Bessel transforms of / are related 
through 



fim{k) ■■ 



If / is band-limited in terms of its Fourier-Laguerre decomposition, the latter 
sum is finite and both transforms can be calculated exactly since jip(k) admit 
the exact analytic formula Eqn. {I2| . 
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on the radial half-line leads to a sampling theorem on B^. 
For a signal with angular and radial band-limits L and P, 
respectively, all of the information content of the signal 
is captured in N = P[(2L - 1)(L - 1) + 1] - 2PL'^ 
samples, yielding an exact Fourier-Laguerre transform on . 
The three-dimensional sampling consists of spherical shells, 
discretised according to a sampling theorem (where here we 
adopt the MW sampling theorem), located at the nodes of 
the radial sampling. The radial sampling may furthermore be 
rescaled to any spherical region of interest [0, i?] x S'^ using 
the parameter r to dilate or contract the radial quadrature rule. 

III. Wavelets on the Ball 

The exactness of the Fourier-Laguerre transform supports 
the design of an exact wavelet transform on the ball. In 
this section we first define a three-dimensional convolution 
operator on the ball, derived from the convolutions defined 
on the sphere and on the radial half-line. We then construct 
flaglets through an exact tiling of Fourier-Laguerre space, 
leading to wavelet kernels which are spatially localised and 
form a tight frame. Furthermore, each kernel projects onto an 
angular scale which is invariant under radial translation. We 
finally introduce a multiresolution algorithm to compute the 
flaglet transform and capture the information of each wavelet 
scale in the minimal number of samples on the ball, while 
optimising the computational cost of the transform. 

A. Convolutions 

The convolution of two functions / and in a (Hilbert) 
space of interest is often defined by the inner product of / 
with a transformed version of h. In standard Fourier analysis 
this transformation is the natural translation. Likewise, for 
two signals on the sphere f^h G the convolution is 

constructed from the rotation operator IZ^j'. 



(21) 



where, here and henceforth, we restrict ourselves to axisym- 
metric kernels h, so that the rotation is only parameterised 
by an angle u = (6>, ^) |[8|, fST]] . The spherical harmonic 
decomposition of fi^h is given by the product of the individual 
transforms: 



{f^h),^ = {f^h\Y,rr^) 



47r 
2ZT1 



fimh^Q, (22) 



with fim = {f\yim) and hmSmO = (^iX^m). 

Similarly, we introduce a translation operator % to con- 
struct the convolution of two functions on the radial half-line 

f,heLHR+y. 

{f*g){r) ^ iflTrh) = [ dr'r'^fir') (Trh)* (r'). (23) 



The convolution in Laguerre space p8|-||6Q| is defined such 
that the action % on the basis functions is 



in which case fi^h simplifies to a product in spherical Laguerre 
space, yielding 



if^h)^ = {f^h\Kp) = fph;, 



(25) 



where fp = {f\Kp) and hp = {h\Kp). Consequently any 
function / which is translated by a distance r on the radial 
half-line has each coefficient fp transformed into fpKp{r). 
This operation corresponds to a translation with a damping 
factor, which is illustrated on a wavelet ker nel in Figure |3] 
(the wavelet kernel itself is defined in Section III-C )^ 

Finally, we define the convolution of two functions on the 
ball fjhe L'^{B^), where h is again assumed to be axisym- 
metric in the angular direction, by combining the convolution 
operators defined on the sphere and radial half-line, yielding 



{f^h){r) = iflTrlZ^h) (26) 
= / dV/(rO(r.7^^/l)*(rO. (27) 
The convolution is given in harmonic space by the product 

(/^^)£mp = {f^h\Zemp) = 



47r ^ 

with fimp = {f\Zemp) and heopSmO = {h\Zimp)- 



(28) 



B. Exact flaglet transform 

With an exact harmonic transform and a convolution oper- 
ator defined on the ball in hand, we are now in a position to 
construct the exact flaglet transform on the ball. For a function 
of interest / G L^(P^), we define its j/-th wavelet coefficient 
W^''' e L^{B^) as the convolution of / with the flaglet (i.e. 
wavelet kernel) G L'^{B^): 

iy^"'(r) = (/^^^■^■')(r) = (/|r.7^^^^■^■'). (29) 

The scales j and / respectively relate to angular and radial 
spaces. Since we restrict ourselves to axisymmetric kernels, 
the wavelet coefficients are given in Fourier-Laguerre space 
by the product 



imp 



Air ■ ■' 

2i '^^'^P ^^p ' 



(30) 



where W^^^ = 



{W"^'' \Zimp), fimp {f\Zimp) and 
^^Qp^mQ — (^•^•^ \Zimp)' The wavelet coefficients contain 
the detail information of the signal only; a scaling function 
and corresponding scaling coefficients must be introduced to 
represent the low-frequency, approximate information of the 
signal. The scaling coefficients G L'^{B^) are defined by 
the convolution of / with the scaling function ^ G L'^{B^): 



l^*(r) = (/*$)(r) = (/|r.7e^$), 



or in Fourier-Laguerre space, 



47r ^ 
2£+ I f^^P^^Opj 



(31) 



(32) 



{TrKp){r')^K;{r)Kp{r'), 



(24) 



where Wg^^ = {W^\Zi^p) and ^mpS^o = {*|^£mp). 

^Note that this translation operator may also be viewed in real space as a 
convolution with a delta function, similarly to the Euclidian convolution. 
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0.1 0.2 0.3 0.4 



(a) Wavelet kernel translated by r = 0.2 





0.1 0.2 0.3 0.4 



(b) Wavelet kernel translated by r = 0.3 





0.1 0.2 0.3 0.4 



(c) Wavelet kernel translated by r = 0.4 

Fig. 3. Slices of an axisymmetric flaglet wavelet kernel constructed on the 
ball of radius R = 1, translated along the radial half-line. The chosen kernel 
has j = j' = 5 and is constructed at resolution P = L = 64. For clarity we 
zoomed on the range r G [0, 0.5] (the sHce hence relates to a ball of radius 
r = 0.5). The three-dimensional wavelet can be visualised by rotating this 
slice around the vertical axis passing through the origin. The translation on 
the radial half-line not only translates the main feature (the wavelet peak) but 
also accounts for a damping factor. Flaglets are well localised in both real and 
Fourier-Laguerre spaces and their angular aperture is invariant under radial 
translation. 



Provided the flaglets and scaling function satisfy an admis- 
sibility property, a function / may be reconstructed exactly 
from its wavelet and scaling coefficients by 



f{r) = I dVH^*(r')(r.7^^$)(r') 



(33) 



J J' 

E E 

j = JO j' = J'Q 



or equivalently in harmonic space by 



47r 

2ZTr 



47r 



J J' 



(34) 



21- 



j = Joj' = Jl^ 



The parameters Jo, Jq, J and J' defining the minimum and 
maximum scales must be defined consistently to extract and 
reconstruct all the information contained in /. They depend 
on the construction of the flaglets and scaling function and are 
defined explicitly in the next section. 

Finally, the admissibility condition under which a band- 
limited function / can be decomposed and reconstructed 
exactly is given by the following resolution of the identity: 



47r 
2ZT1 



J' 



- E E \< 



= 1, 



(35) 



We may now construct flaglets and scaling functions that 
satisfy this admissibility property and thus lead to an exact 
wavelet transform on the ball. 

C. Flaglets and scaling functions 

We extend the notion of harmonic tiling |[8|, ||9|, | [6T| 
to the Fourier-Laguerre space and construct axisymmetric 
wavelets (flaglets) well localised in both real and Fourier- 
Laguerre spaces. We first define the flaglet and scaling function 
generating functions, before defining the flaglets and scaling 
function themselves. 

We start by considering the Schwartz function with 
compact support 



s{t) 



e 

0, 



t e 



-1,1] 
-1,1] 



(36) 



for t G M. We introduce the positive real parameter A G 
to map s(t) to 

2A 



A- 1 



{t - 1/A) - 1 



(37) 



which has compact support in [;^,1]- We then define the 
smoothly decreasing function k\ by 



kx{t) 



(38) 



which is unity for t < 1/A, zero for t > 1, and is 
smoothly decreasing from unity to zero for t G [1/A, 1]. 
Axisymmetric flaglets are constructed in a two-dimensional 
space corresponding to the harmonic indices i and p. We 
associate A with ^- space and we introduce a second parameter 
u associated with p- space, with the corresponding functions 
Sj^ and kj^. We define the flaglet generating function by 



i^x{t) = y^kx{t/X)-kx{t) 
and the scaling function generating function by 



(39) 



(40) 



with similar expressions for and ?7^, complemented with a 
hybrid scaling function generating function 



kx{t/X)K{t') 
■ kx{t)K{t'/v) 
kx{t)K{t') 



(41) 



1/2 
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The flaglets and scaling function are constructed from their 
generating functions to satisfy the admissibihty condition 
given by Eqn. jisj). A natural approach is to define "^^i^p 
from the generating functions hix and k,j^ to have support on 
[Ai-i,A^'+i] X yielding 



^ imp 



2£ + l 
47r 



(42) 



With these kernels, Eqn. (35) is satisfied for i > and 
p > w^o, where Jq and Jq are the lowest wavelet scales used 
in the decomposition. The scaling function <l> is constructed 
to extract the modes that cannot be probed by the flaglets 



elsewhere. 



To satisfy exact reconstruction, J and J' are defined from the 
band-limits by J = [log;,(L - 1)] and J' = \\og^{P - 1)1. 
The choice of Jo and Jq is arbitrary, provided that < Jo < J 
and < Jq < J^ This framework generalises the notion of 
the harmonic tiling used to construct exact wavelets on the 
sphere Q, ||9j; in fact, the flaglets defined here reduce in 
angular part to the wavelets defined in 1 8 1 for the axisymmetric 
case. The flaglets and scaling function tiling of the Fourier- 
Laguerre space of the ball is illustrated in Figure |4] Flaglets 
and the scaling function may be reconstructed in the spatial 
domain from their harmonic coefficients. In Figure [5] flaglets 
are plotted in the spatial domain for a range of different scales; 
translated flaglets are plotted in Figure [3] The flaglets are well 
localised in both real and Fourier-Laguerre spaces and their 
angular aperture is invariant under radial translation. 

IV. MULTIRESOLUTION ALGORITHM 

In this section we discuss our implementation of the Fourier- 
Laguerre and flaglet transforms. We notably introduce a mul- 
tiresolution algorithm for the flaglet transform to capture each 
wavelet scale in the minimal number of samples on the ball, 
thereby reducing the computational cost of the transform. 
We finally provide accuracy and complexity tests for our 
implementation of both transforms, which we make publicly 
available. 

A. Algorithm 

In our framework, each flaglet ^-^-^ has compact sup- 
port in Fourier-Laguerre space on £ x p e [A-^~^, A-^+^] x 
[u^ ~^,z^^ +^], as shown in Figure [4] Thus, ^^-^ has band- 
limits in £ and p of A-^+^ and +^ respectively. For a band- 
limited function / G L'^{B^), recall that the j/-th wavelet 



contribution is given by the simple product of Eqn. (30) in 
harmonic space. Consequently, the band-limits of W^^^ are 
given by the minimum of the band-limits of / and . 

^Note that despite its piecewise definition ^£rnp is continuous along and 
across the boundaries p = ly^o and i = A'^o . 




Fig. 4. Tiling of Fourier-Laguerre space at resolution L = = 64 for flaglet 
parameters X = u = 2, giving J = J' = 7 . Flaglets divide Fourier-Laguerre 
space into regions corresponding to specific scales in angular and radial space. 
The scaling part, here chosen as Jq = Jq = 4, is introduced to cover the 
low frequency region and insures that large scales are also represented by the 
transform. 




(a) (i,i') = (4,5) 



(b) (i,i') = (4,6) 





(c) (i,i') = (5,5) 



(d) (i,/) = (5,6) 



Fig. 5. Slices of four successive axisymmetric flaglet wavelet kernels, probing 
different scales in angular and radial space. The flaglet parameters are A = 
u = 2 and the kernels are constructed at resolution L = A^ = 92ona ball of 
radius R = 1. For visualisation purposes we show the flaglets corresponding 
to j G {4, 5} and j' G {5,6}, translated to r = 0.3 and zoomed on the range 
r G [0, 0.4]. Kernels of angular order j = 4 (first row) probe large angular 
scales compared to those of order j = 5 (second row). Similarly, kernels of 
radial order j' = 5 (first column) probe large radial scales compared to those 
of order j' = 5 (second column). 
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Thus, for j < J ox j' < J' the wavelet scale W^^^ 
can be represented in fewer samples than /, without any 
loss of information. We exploit this property by designing 
a multiresolution approach where each wavelet scale is rep- 
resented in real space with the smallest number of samples 
necessary. Note that the scaling function must be used at full 
resolution since its angular and radial band-limits are L and 
P respectively. To summarise the multiresolution algorithm, 
although / is decomposed at full resolution, the wavelets 
coefficients are reconstructed in real space with the minimum 
number of samples supporting their band-limits. This leads 
to a significant reduction in computation time, which is then 
dominated by the small number of full resolution Fourier- 
Laguerre transforms. 

B. Fast implementation 

Our implementation of the algorithms of this article is made 
available in the following three packages, which are written in 
C and include MAT LAB interfaces for most high-level features, 
and are described in turn: 

• FLAG: spherical Laguerre transform and Fourier- 
Laguerre transforms on the ball (exact spherical Bessel 
and Fourier-Bessel decompositions are optional features 
that additionally require the GNU Math Librar}Q. 

• S2LET: axisymmetric wavelet transform on the sphere 
through harmonic tiling. 

• FL AGLET: axisymmetric flaglet transform on the ball, 
combining FLAG and S2LET to construct flaglets in 
Fourier-Laguerre space through harmonic tiling. 

We make these three packages publicly available]^ All pack- 
ages require SSHt|^ which implements fast and exact algo- 
rithms to perform the forward and inverse spherical harmonic 
transforms correspondin g to the MW sampling theorem |56|. 
SSHT requires the FFTv|^ package. 

Since the naive spherical harmonic transform scales as 
0{L^) and the spherical Laguerre transform scales as 0{P'^), 
the naive complexity of the Fourier-Laguerre transform 
is 0{P'^L^). However, rather than computing triple inte- 
grals/sums over the ball directly, it is straightforward to 
show that the Fourier-Laguerre transform can be performed 
separately on the sphere and on the radial half-line, like the 
Fourier-Bessel transform |50|. Since the angular and radial 
samplings are separable, the related transforms can be com- 
puted independently through a separation of variables, so that 
the complexity reduces to 0{Q^) for Q ~ P ~ L. The 
separation of variables also means we are able to exploit high- 
performance recurrences and algorithms that exist for both 
the spherical Laguerre and spherical harmonic transforms. In 
particular, the radial basis functions Kp{r) are calculated using 
a normalised recurrence formula derived from the recurrence 
on the Laguerre polynomials. Moreover, a critical point for the 
accuracy of the Fourier-Laguerre transform is the computation 

^http://www.gnu.org/software/gsl/ 

\http://www.flaglets.org/ 

%ttp://w ww.j asonmce wen.org/| 
ihttp://www. ff tw. org/| 



of the Gauss-Laguerre quadrature, for which we use the previ- 
ous normalised recurrence complemented with an appropriate 
root-finder algorithm. The fast spherical harmonic transforms 
implemented in the SSHT package use the Trapani & Navaza 
method |62| to efficiently compute Wigner functions (which 
are closely related to the spherical harmonics) through re- 
cursion]^ These fast spherical harmonic transform algorithms 
1 56 1 scale as 0{L^). The final complexity achieved by the 
Fourier-Laguerre transform is thus 0{Q^). 

The flaglet transform (forward and inverse) is calculated in 
a straightforward manner in Fourier-Laguerre space, thus its 
computation is dominated by the Fourier-Laguerre transform 
of the signal, approximation coefficients, and wavelets coeffi- 
cients at all scales, requiring [( J + 1 — Jq){J' + 1 — Jq) + 2] 
Fourier-Laguerre transforms. If all wavelet contributions are 
reconstructed at full resolution in real space, the overall 
wavelet transform scales as 0([{J + 1 — Jq){J' + 1 — 
Jo)+2]Q^). Note that J and J' depend on the band-limits L 
and P and the parameters A and respectively. However, in 
the previous section we established a multiresolution algorithm 
that takes advantage of the band-limits of the individual 
flaglets. With this algorithm, only the scaling function and the 
finest wavelet scales (i.e. j G { J — 1, J} and / G { — 1, J'}) 
are computed at maximal resolution corresponding to band- 
limits L and P. The complexity of the overall multiresolution 
flaglet transform is then dominated by these operations and 
scales as 0{Q^). 

C. Numerical validation 

In this section we evaluate FLAG and FLAGLET in terms 
of accuracy and complexity. We show that they achieve 
floating-point precision and scale as detailed in the previous 
section. In both cases we consider band-limits L = P = 2* 
with i e {2, . . . , 9} and generate sets of harmonic coef- 
ficients firnp following independent Gaussian distributions 
A/'(0, 1). We then perform either the Fourier-Laguerre or 
the flaglet decomposition, before reconstructing the har- 
monic coefficients, therefore denoted by f/^p. We evalu- 
ate the accuracy of the transforms using the error metric 
e = max | fi^np — ff^p I ' which is theoretically zero for both 
transforms since all signals are band-limited by construction. 
The complexity is quantified by observing how the compu- 
tation time tc = [Synthesis + Analysis] /2 scalcs with the band- 
limits, where the synthesis and analysis computation times, 
^synthesis ^nd Analysis respectively, are defined explicitly for the 
two transforms in the paragraphs that follow. The stability 
of both e and is checked by averaging over hundreds of 
realisations of /^^^ in the cases i G {2, . . . , 7} and a small 
number of realisations for i G {8,9}. Recall that for given 
band-limits L and P the number of samples on the ball 
required by the exact quadrature isA/' = P[(2L— 1)(L— 1)+1]. 
All tests were run on a 2.5GHz Core 15 processor with 8GB 
of RAM. 

The results of these tests for the Fourier-Laguerre transform 
are presented on Figure [6] The indicators e and tc are plotted 

^^Alternatively, Risbo's method could also be used to compute Wigner 
functions (631. 
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against the number of samples N. Each test starts from 
coefficients f£mp randomly generated. The synthesis refers to 
constructing the band-limited signal / from the decomposition 
fimp- The analysis then corresponds to decomposing / into 
Fourier-Laguerre coefficients f/^p. As shown in Figure [61 
FLAG achieves very good numerical accuracy, with numerical 
errors comparable to floating-point precision, and computation 
time scales as 0{Q^), in agreement with theory. 

The results of similar tests for the flaglet transform (entirely 
performed in real space) are presented on Figure [7] As 
previously, the indicators e and tc are plotted against the 
number of samples N. Since we evaluate the flaglet transform 
in real space, a preliminary step is required to construct a 
band-limited signal / from the randomly generated fimp- 
This step is not included in the computation time since its 
only purpose is to generate a valid band-limited test signal 
in real space. The analysis then denotes the decomposition 
of / into wavelet coefficients W^^^ and scaling coefficients 
on the ball. The synthesis refers to recovering the signal 
yrec f^Q^ thesc Coefficients. The final step, which is not 
included in the computation time, is to decompose f^^ into 
Fourier-Laguerre coefficients /^^^ in order to compare them 
with fimp- As shown in Figure ItI FLAGLET achieves very 
good numerical accuracy, with numerical errors comparable 
to floating-point precision. Moreover, the full resolution and 
multiresolution algorithms are indistinguishable in terms of 
accuracy. However, the latter is ten times faster than the 
former since only the scaling function and a small number 
of wavelet coefficients are computed at full resolution. As 
shown in Figure [t] computation time scales as 0{Q^) for 
both algorithms, in agreement with theory. 

V. Denoising Illustration 

In this section we illustrate the use of the flaglet transform in 
the context of a simple denoising problem. We consider two 
datasets naturally defined on the ball and contaminate them 
with band-limited noise. We compute the flaglet transform 
of the noisy signal and perform simple denoising by hard- 
threshold the wavelet coefficients. We reconstruct the signal 
from the thresholded wavelet coefficients and examine the 
improvement in signal fidelity. 

A. Wavelet denoising 

Consider the noisy signal y = s -\- n G L'^{B^), where 
the signal of interest s G L'^{B^) is contaminated with noise 
n G L'^{B^). A simple way to evaluate the fidelity of the 
observed signal y is to examine the signal-to-noise ratio, which 
we define on the ball by 

SNR(y) = 101ogio ||J'^"L - (43) 

\\y ^\\2 

The signal energy is given by 




where the final equality follows from a Parseval relation on 
the ball (which follows directly from the orthogonality of the 




2? 2'^^ 2'^^ 2^^ 2^^ 2^^ 2^^ 

N 

(a) Numerical accuracy of the Fourier-Laguerre transform 



10' 




2? 2^'^ 2^*^ 2P 2^^ 2^^ 2^^ 

(b) Computation time of the Fourier-Laguerre transform 



Fig. 6. Numerical accuracy and computation time of the Fourier-Laguerre 
transform computed with FLAG, where corresponds to the number of 
samples on the ball required to capture all the information contained in the 
band-limited test signal. We consider L = P = 2^ with i G {2,. ..,9}. 
These results are averaged over many realisations of random band-limited 
signals and were found to be very stable. Very good numerical accuracy is 
achieved, with numerical errors comparable to floating-point precision, found 
empirically to scale as O(Q^) as shown by the red line in panel (a), where 
Q P L. Computation time scales as O(Q^) as shown by the red line 
in panel (b), in agreement with theory. 

Fourier-Laguerre basis functions). In practice, we compute 
signal energies through the final Fourier-Laguerre space ex- 
pression to avoid the necessity of an explicit quadrature rule. 

We seek a denoised version of y, denoted hy d e L'^{B^), 
such that SNR((i) is as large as possible in order to extract the 
informative signal s. We take the flaglet transform of the noisy 
signal since we intend to denoise the signal in wavelet space, 
where we expect the energy of the informative signal to be 
concentrated in a small number of wavelet coefficients while 
the noise energy will be spread over many wavelet coefficients. 
Since the flaglet transform is linear, the wavelet coefficients 
of the j/-th scale of the noisy signal is simply the sum of the 
individual contributions: 

yn (^) = gjf (^) ^ j^jf (^)^ (45) 

where capital letters denote the wavelet coefficients, i.e. 
yn' =y^ =si, y^n' and N^^' = ni. 

In the illustrations performed here, we assume the noise 
model 

E(l^^mpP) = [pj ^U'^mm'Spp', (46) 

which corresponds to a white noise for the angular space 
with a dependence on the radial mode p, where E(-) denotes 
ensemble averages. We do not opt for a white noise in radial 
space (i.e. E(|n£^pp) = cr'^Sii'Smm'^pp') because the latter 
has its energy concentrated in the centre of the ball due 
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(a) Numerical accuracy of the flaglet transform 




(b) Computation time of the flaglet transform 

Fig. 7. Numerical accuracy and computation time of the flaglet transform 
computed with FLAGLET, where corresponds to the number of samples 
on the ball required to capture all the information contained in the band- 
limited test signal. We consider L = P = 2^ with i G {2,... ,9}, with 
parameters X = u = 2, Jq = Jq = 0. These results are averaged over 
many realisations of random band-limited signals and were found to be very 
stable. The flaglet transform is either performed at full-resolution (dashed 
lines) or with the multiresolution algorithm (solid lines). Very good numerical 
accuracy is achieved by both the full resolution and multiresoltion algorithms 
(which achieve indistinguishable accuracy), with numerical errors comparable 
to floating-point precision, found empirically to scale as 0(PL) as shown by 
the red line in panel (a). The multiresolution algorithm is ten times faster than 
the full-resolution approach. Computation time scales as O(PL^) for both 
algorithms as shown by the red line in panel (b), in agreement with theory. 



to the shape of the spherical Laguerre basis functions. The 
p-dependence gives a greater weight to small-scale radial 
features and hence yields a more homogeneous noise on the 
ball, which is more useful for visualisation purposes. For this 
noise model one can show that the expected covariance of the 
wavelet coefficients of the j/-th scale reads 



Denoising is performed by hard-thresholding the 
wavelet coefficients Y^^ , where the threshold is taken 
as T(r) = 3a^^ (r). The wavelet coefficients of the denoised 
signal D^^ = di^ are thus given by 



0, if Y^^ (r) < T(r) 

Y^^ (r), otherwise. 



(48) 



The denoised signal d is then reconstructed from its wavelet 
coefficients and scaling coefficients (the latter are not thresh- 
olded and thus not altered). To assess the effectiveness of this 
simple flaglet denoising strategy when the informative signal 
s is known, we compute the SNR of the denoised signal 
and compare it to the SNR of the original noisy signal. In 



what follows we apply this simple denoising technique to two 
datasets naturally defined on the ball. 

B. Examples 

The first dataset we consider is the full-sky Horizon sim- 
ulation |64| : an N-body simulation covering a IGpc periodic 
box of 70 billion dark matter particles generated from the 
concordance model cosmology derived from 3 -year Wilkinson 
Microwave Anisotropy Probe (WMAP) observations | [65| . The 
purpose of such a simulation is to reproduce the action of 
gravity (and to a minor extent galaxy formation) on a large 
system of particles, with the initial conditions drawn from a 
cosmological model of interest. The outcome is commonly 
used to confront astrophysical models with observations. For 
simplicity we only consider a ball of IMPc radius centered at 
the origin so that the structures are of reasonable size. Figure [5] 
shows the initial data, band-limited at L = P = 128, as well 
as their wavelet coefficients with \ = v = 2, J = J' = 7 
and scaling coefficients for Jq = = 6 since the lower scale 
indices do not contain a great deal of information. We see 
that the filamentary distribution of matter is naturally suited 
to a flaglet analysis on the ball since the informative signal 
is likely to be contained in a reduced number of wavelet 
coefficients. The original data are corrupted by the addition 



of random noise defined by Eqn. (46) for an SNR of 5dB. 
The wavelet denoising procedure described previously is then 
applied. The denoised signal is recovered with an SNR of 
lldB, highlighting the effectiveness of this very simple flaglet 
denoising strategy on the ball. The results of this denoising 
illustration are presented in Figure [9] 

The second dataset we consider is Ritsema's seismological 
Earth model of shear wavespeed perturbations in the mantle, 
known as S40RTS (13), |14|, \66\^ The model suppHes 
spherical harmonic coefficients in the angular dimension and 
radial spline coefficients in the depth dimension to define a 
signal on the ball, which we band-limit. Contrarily to the 
first example, Ritsema's model does not contain a lot of 
structure at the smallest scales but essentially contains large- 
scale features. As previously, the original data are corrupted by 



the addition of random noise defined by Eqn. (46) for an SNR 
of 5dB. The flaglet denoising procedure described previously 
is then applied. The denoised signal is recovered with an 
SNR of 17dB, again highlighting the effectiveness of this very 
simple flaglet denoising strategy on the ball. As expected, the 
improvement in SNR is better than for the previous dataset 
since the informative signal is mainly captured by a few large 
wavelet scales. The results of this denoising illustration are 



presented in Figure 10 



VI. Conclusions 

One reason an exact wavelet transform of a band-limited 
signal on the ball has not yet been derived is due to the absence 
of an exact harmonic transform on the ball. We have taken 
advantage of the orthogonality of the Laguerre polynomials on 
R+ to define the spherical Laguerre transform, a novel radial 



^ ^tp://www.earth.lsa.umich.edu/"^jritsema/ 
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(a) Band-limited data (b) Scaling coefficients (a) Band-limited data (b) Noise 




(c) 0',/) = (6,6) (d) a,/) = (6,7) 




(e) 0',/) = (7,6) (f) (jJ') = (7J) 



Fig. 8. Flaglet decomposition of the N-body simulation dataset considered 
for the first denoising example. The initial dataset was pixelised and band- 
limited ai L = P = 12S. The flaglet parameters arc X = u = 2 (giving 
J = J' = 7) and the scahng coefficients correspond to Jo = Jq = 6 
since the lower scale indices do not contain a great deal of information. The 
four wavelet coefficients together with the scaling coefficients decompose 
the initial dataset exactly, i.e. the original signal can be recovered perfectly 
from these wavelet and scahng coefficients. All signals were oversampled on 
L = P = 256 for visualisation purposes. 

transform that admits an exact quadrature rule. Combined with 
the spherical harmonics, we used this to derive a sampling 
theorem and an exact harmonic transform on the ball, which 
we call the Fourier-Laguerre transform. A function that is 
band-limited in Fourier-Laguerre space can be decomposed 
and reconstructed at floating-point precision, and its Fourier- 
Bessel transform can be calculated exactly. For radial and an- 
gular band-limits P and L, respectively, the sampling theorem 
guarantees that all the information of the band-limited signal 
is captured in a finite set of N = P[{2L - 1)(L - 1) + 1] 
samples on the ball. 

We have developed an exact wavelet transform on the ball, 
the so-called flaglet transform, through a tiling of the Fourier- 
Laguerre space. The resulting flaglets form a tight frame and 




(c) Noisy signal (d) Denoised signal 



Fig. 9. Denoising of an N-body simulation. The data are contaminated with 
a band- limited noise and decomposed into wavelet coefficients. Denoising is 
performed by a simple hard-thresholding of the wavelet coefficients, following 
a noise model. The denoised signal is then reconstructed from the thresholded 
wavelet coefficients. In this example, for an initial SNR of 5dB, the flaglet 
denoised signal is recovered with SNR of SNR = lldB (with resolution 
L = P = 128, oversampled on L = P = 256 and using flaglet parameters 
X = u = 2, Jo = Jl^ = 0, giving J = J' = 7). 



are well localised in both real and Fourier-Laguerre spaces. 
Their angular aperture is invariant under radial translation. 
We furthermore established a multiresolution algorithm to 
compute the flaglet transform, capturing all the information 
contained in each wavelet scale in the minimal number of 
samples on the ball, thereby reducing the computation cost of 
the flaglet transform considerably. 

Flaglets are a promising new tool for analysing signals on 
the ball, particularly for extracting spatially localised features 
at different scales of interest. Exactness of both the Fourier- 
Laguerre and the flaglet transforms guarantees that any band- 
limited signal can be analysed and decomposed into wavelet 
coefficients and then reconstructed without any loss of in- 
formation. To illustrate these capabilities, we considered the 
denoising of two different datasets which were contaminated 
with synthetic noise. A very simple flaglet denoising strategy 
was performed by hard-thresholding the wavelet coefficients 
of the noisy signal, before reconstructing the denoised signal 
from the thresholded wavelet coefficients. In these illustrations 
a considerable improvement in SNR was realised by this sim- 
ple flaglet denoising strategy, demonstrating the effectiveness 
of flaglets for the analysis of data defined on the ball. Our 
implementation of all of the transforms and examples detailed 
in this article is made publicly available. In future work we 
intend to revoke the axisymmetric constraint by developing 
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(a) Band-limited data (b) Noise 




(c) Noisy signal (d) Denoised signal 



Fig. 10. Denoising of a seismological Earth model. The data are contaminated 
with a band-limited noise and decomposed into wavelet coefficients. Denois- 
ing is performed by a simple hard-thresholding of the wavelet coefficients, 
following a noise model. The denoised signal is then reconstructed from the 
thresholded wavelet coefficients. In this example, for an initial SNR of 5dB, 
the flaglet denoised signal is recovered with SNR of 17dB (with resolution 
L = P = 128 and using flaglet parameters X = v = 3, Jo = Jq = 0, 
giving J = J' = 7). 



directional flaglets. 
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